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ABSTRACT 

We present the results of a systematic survey of numerical solutions to the 
coagulation equation for a rate coefficient of the form A\j oc {i tl j v + i v J M ) and 
monodisperse initial conditions. The results confirm that there are three classes of 
rate coefficients with qualitatively different solutions. For v < 1 and A = \i-\-v < 1, 
the numerical solution evolves in an orderly fashion and tends toward a self-similar 
solution at large time t. The properties of the numerical solution in the scaling 
limit agree with the analytic predictions of van Dongen and Ernst. In particular, 
for the subset with \i > and A < 1, we disagree with Krivitsky and find that 
the scaling function approaches the analytically predicted power-law behavior at 
small mass, but in a damped oscillatory fashion that was not known previously. 
For v < 1 and A > 1, the numerical solution tends toward a self-similar solution 
as t approaches a finite time £q. The mass spectrum n& develops at £q a power- 
law tail nk oc k~ T at large mass that violates mass conservation, and runaway 
growth/gelation is expected to start at £ cr it = to in the limit the initial number of 
particles no — > oo. The exponent r is in general less than the analytic prediction 
(A + 3)/2, and t = K/[(X - l)n A n ] with K = 1-2 if A > 1.1. For v > 1, the 
behaviors of the numerical solution are similar to those found in a previous paper 
by us. They strongly suggest that there are no self-consistent solutions at any time 
and that runaway growth is instantaneous in the limit no — > oo. They also indicate 
that the time t cr it for the onset of runaway growth decreases slowly toward zero 
with increasing no- 
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1. INTRODUCTION 



Smoluchowski's coagulation equation is the mean-rate equation that describes 
the evolution of the mass spectrum of a collection of particles due to successive 
mergers. It is widely used for modeling growth in many fields of science. Examples 
include planetesimal accumulation, mergers in dense clusters of stars, coalescence 
of interstellar dust grains, and galaxy mergers in astrophysics, aerosol coalescence 
in atmospheric physics, colloids, and polymerization and gelation (see, e.g., Drake 
1972, Ernst 1986, Jullien and Botet 1987, Lee 1993, 2000, and references therein). 

If the masses of the particles are integral multiples of a minimum mass mi, the 
coagulation equation is written in discrete form as 



where rik is the number of particles of mass = km\ in a volume V and Aij is 
the rate coefficient (or coagulation kernel) for mergers between particles of mass rrii 
and rrij. 1 In equation (1), it is assumed that the merging of two particles of mass 
rrii and rrij results in one particle of mass rrii + rrij. The coagulation equation can 
also be written in continuous form as 



where n(m)dm is the number of particles of mass between m and m + dm and A mjTn i 
is the rate coefficient for mergers between particles of mass m and m' . 

Examples of the rate coefficient A^ as a function of rrii and rrij (or equivalently 
% and j) that arise in various problems can be found in the references cited above. 
Most rate coefficients used in the literature are homogeneous functions of degree 
A, i.e., A aiy aj = a x Aij. The exponent A specifies the mass dependence of the 
probability of merger for two particles of comparable mass (i ~ j). It is also useful 
to classify A^ according to the exponents fi and v for the merger between a light 
particle and a heavy particle: A^ oc i^j v for i <C j and \i + v = A (see, e.g., Ernst 
1986). For example, A^ oc % + j has n = 0, v = 1, and A = 1. 

For a few simple rate coefficients and monodisperse initial conditions (i.e., no 
particles of mass mi at t = 0), there are exact analytic solutions to the discrete form 

1 We can interpret rik as the concentration (i.e., the number of particles per unit 
volume) if we replace A^ by A'- = V A^. 
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of the coagulation equation (Trubnikov 1971, Ernst 1986, and references therein). 
The analytic solutions for A^ oc constant and i + j show orderly evolution of a 
smooth mass spectrum at all times, and they agree with the results from Monte 
Carlo simulations of the merger process in the limit no — > oo (with uq/V fixed). 
These two cases are examples of orderly growth. The analytic solution for A^ oc ij 
develops a power-law tail with nk oc /c -5 / 2 at large k as t | to = l/(no^4n). This 
power-law tail violates mass conservation because it implies a nonzero mass flux to 
the infinite-mass bin. In this case, the results from Monte Carlo simulations in the 
limit fi-o — ► oo agree with the solution to the coagulation equation at t < to, but 
they show a transition from a smooth mass spectrum to a smooth spectrum plus 
a massive runaway particle at t = to (Spouge 1985, Wetherill 1990). The runaway 
particle (gel) acquires a mass much larger than that of the other particles (sol) 
in the system and becomes detached from the smooth mass spectrum of the rest 
of the particles at t > to- This phenomenon is known as runaway growth in the 
astrophysics literature, and the transition is considered to be the gelation transition 
in studies of polymerization and gelation. 

For most rate coefficients, there are no exact analytic solutions to the 
coagulation equation. However, there are extensive analytic results on the 
asymptotic properties of the solutions for A^ with v < 1 (see, e.g., review by Ernst 
1986). It is important to note that some of these analytic results (such as the shape 
of the mass spectrum at small and large mass) are derived based on assumptions 
(such as self-similar evolution) that have not been verified. Nevertheless, the 
analytic results indicate that there are qualitatively two types of solutions to the 
coagulation equation if v < 1: 

(1) if v < 1 and A < 1, the solution shows orderly growth at all times; 

(2) if v < 1 and A > 1, the solution develops in a finite time to a power-law tail at 
large mass that violates mass conservation. 

For A^ with essentially v < 1, it has been proved that a solution to the coagulation 
equation exists for all times (including t > to if A > 1) and that the coagulation 
equation is the limit of finite system (whether or not the runaway particle and the 
other particles are allowed to interact at t > to if A > 1) (Leyvraz and Tschudi 
1981, Spouge 1985, Bak and Heilmann 1994, Jeon 1998). (For v = 1 and A > 1, we 
have the example A^ oc ij, where the coagulation equation needs to be modified 
for t > to if there is sol-gel interaction; Ziff et al 1983, Bak and Heilmann 1994.) 

Several authors have investigated the properties of the solutions to the 
coagulation equation for Aij with v > 1, using series expansion of the mass spectrum 
nk{t) about t = and moments of the mass spectrum (McLeod 1962, Hendriks et 
al 1983, Ernst et al 1984, van Dongen 1987a). The results suggest that 

(3) if v > 1, there are no self-consistent solutions that conserve mass at any time. 
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An alternative to the analytic approach is to solve the coagulation equation 
numerically. In Paper I (Lee 2000), a numerical code that can yield accurate 
solutions to the discrete form of the coagulation equation, equation (1), with a 
reasonable number of numerical mass bins was developed. The numerical code was 
used to study solutions to the coagulation equation for that are limiting cases 
for gravitational interaction. We considered geometric or gravitational focusing 
dominated cross-section, mass-independent or equipartition velocity dispersion, and 
the power-law index of the mass-radius relation /3 = 1/3 (for planetesimals) or 2/3 
(for stars). For the two cases with geometric cross-section and (3 = 1/3, which 
have v < 1 and A < 1, the mass spectrum evolves in an orderly fashion and tends 
towards a self-similar solution at large time. For the remaining cases, which have 
v > 1, the numerical mass spectrum shows, after some evolution, an exponential 
drop in an intermediate mass range and a power-law tail of the form njt oc k~ u 
(or n(m) oc m~ v ) at large mass. This mass spectrum is not self-consistent because 
the power-law tail implies a mass flux 2 and, if 1 < v < 2, a cumulative mass that 
diverge with the maximum particle mass, m max , included in the computational 
grid. The time at which the power-law tail develops decreases toward zero as the 
numerical parameter n m i n decreases (see Section 3 for the definition of n m i n ). Thus 
the numerical results strongly suggest that there are no self-consistent solutions to 
the coagulation equation at any time if v > 1. We also considered a case with (3 = 
as an example with v < 1 and A > 1, and its mass spectrum develops a power-law 
tail that violates mass conservation in a finite time t . We discussed a simplified 
merger problem that illustrates the qualitative differences in the solutions to the 
coagulation equation for the three classes of Aij. The results in Paper I (and the 
analytic results cited above) strongly suggest that there are two types of runaway 
growth. For Aij with v < 1 and A > 1, runaway growth starts at a finite time 
^crit = ^o? the time at which the coagulation equation solution begins to violate 
mass conservation, in the limit no — > oo. For A^ with v > 1, runaway growth is 
instantaneous in the limit no — >• oo, and there are indications (since decreasing n m ; n 

2 As we pointed out in Paper I, in these cases, a massive particle grows mainly 
by accumulating low-mass particles because of the much larger number of low-mass 
particles. So the growth rate of a massive particle is m = J dm' ' A m ^ m tn{m'}m! oc 
m v , since the integral is dominated by the range m' <C m. Hence nm at the high- 
mass end of the mass spectrum is non-zero and independent of m if n oc m~ v . 
However, we did not point out that the mass flux from particles of mass m' < m to 
particles of mass m' > m is F m ss mn(m)m(m) + J" mmax dm'n(m')m(m'). With nm 
independent of m, F m m nmm m ^ which is independent of m but increases with 
m max - We have verified this by an explicit evaluation of the mass flux (equation 
[12]) for the numerical solutions. 
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is similar to increasing no) that the time t cr it, in units of l/(no^4n), for the onset 
of runaway growth decreases slowly toward zero with increasing initial number of 
particles no- Recent Monte Carlo simulations have shown that the time for all 
particles to coalesce into a single particle decreases as a power of the logarithm of 
no if v > 1 (Malyshkin and Goodman 2001; see also Spouge 1985, Jeon 1999). 

The study in Paper I was focused on the rate coefficient for gravitational 
interaction, and the range of fx and v studied was limited. In particular, fi was 
limited to and ±1/2, and the region \x > and v < 1 was not studied at all. 
Other authors have obtained numerical solutions to the coagulation equation for 
rate coefficients that arise in specific problems. However, we are not aware of 
any study that has systematically surveyed the properties of numerical solutions 
as a function of [i and v (and A) and compared them to the analytic results on 
the asymptotic properties. For such a study, it is important that the numerical 
code used can follow the evolution of the mass spectrum accurately for a long 
time. After the computation for this paper was nearly complete, it came to our 
attention that Krivitsky (1995) has obtained numerical solutions to the continuous 
form of the coagulation equation, equation (2), for A m m / oc (mm') A / 2 (which has 
/j, = v = A/2) and A m ^ m i oc (m + m') v (which has /i = 0). For A mjm > oc (mm') A / 2 
with A < 1, Krivitsky found that the numerical solutions are self-similar at large 
time, but that unlike the analytic result, the asymptotic behavior of the scaling 
function at small mass is not a power law. As we shall see, the latter result is 
incorrect because Krivitsky did not evolve the numerical solutions for a sufficiently 
long time to see the true asymptotic behavior at small mass. The scaling function 
does in fact approach the analytically predicted power law at small mass, but in 
a damped oscillatory fashion that was not known previously. It is unlikely that 
Krivitsky could have solved the coagulation equation accurately for the necessary 
amount of time because the numerical code used by Krivitsky does not conserve 
mass. For A m ^ m i oc (m + m') 1 ' with v > 1, Krivitsky found that the mass spectrum 
develops a slowly decreasing tail at very small time. Only the numerical solution for 
v = 2 was shown. Its evolution is qualitatively similar to that found in Paper I for 
v > 1, but it is not clear that the tail is power-law in nature because the maximum 
particle mass included in Krivitsky's computational grid is not large enough. It 
was also not demonstrated that the numerical solution is not sensitive to the other 
numerical parameters. 

In this paper we present the results of a systematic survey of numerical solutions 
to the coagulation equation. The purpose of this survey is (1) to confirm that there 
are three classes of rate coefficients with qualitatively different solutions to the 
coagulation equation and that the boundaries of these three classes are as stated 
above; (2) to investigate, in the cases where self-consistent solutions exist, whether 
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the solutions approach self-similar solutions as t — > oo or t ] t and whether the 
scaling behaviors agree with the analytic results; (3) to study the dependence of to 
on the exponents /i, u, and A for the runaway growth cases with v < 1 and A > 1; 
and (4) to investigate whether the behaviors of the numerical solutions found in 
Paper I for the cases with v > 1 are valid in general. In Section 2 we describe 
the rate coefficient and the initial conditions used in this survey. In Section 3 
we provide a brief summary of the numerical methods developed in Paper I for 
solving the coagulation equation and additional information on the accuracy of the 
numerical results. The results are presented in Section 4, and the conclusions are 
summarized in Section 5. 

2. RATE COEFFICIENT AND INITIAL CONDITIONS 

In this paper we consider a rate coefficient of the form 

Ay = I (<T + i v f) (3) 

with /j, < v. Note that A ai ^ a j = a^ +iy Aij and Aij oc i^j" for i <C j, consistent with 
the definitions of the exponents fi, is, and A = \i + v in Section 1. Since the rate 
coefficient in equation (3) contains the exponents \i and v as parameters explicitly, 
we can survey the entire (//, v) space by varying \i and v. This rate coefficient 
includes A tJ = (zj) A/2 (for \i = v = A/2) and A i5 = {i v + j u )/2 (for \i = 0), 
which have been used to model polymerization (e.g., Hendriks et al 1983), and 
the cases (/i, v) = (1,4/3) and (2,2), which have been used to model planetesimal 
accumulation and stellar merger (Malyshkin and Goodman 2001). It also has the 
nice property that it includes the three cases with exact analytic solutions to the 
coagulation equation: Ay = 1, (i + j)/2, and ij for v) = (0,0), (0, 1), and (1, 1), 
respectively. 

We have obtained numerical solutions to the discrete form of the coagulation 
equation for the cases shown in figure 1. The ranges of \i and v considered contain 
most of the values encountered in practical applications (but usually for other 
forms of A^). The numerical solutions were computed for the monodisperse initial 
conditions with no particles of mass mi, i.e., nk(t = 0) = noSki, where Ski is the 
Kronecker delta. Hereafter, we adopt units such that no = 1, mi = 1, and An = 1. 
With this set of units, rrik = k and time is in units of l/(noAn), the timescale for 
every particle of mass mi to merge with another particle of mass mi. 

3. NUMERICAL METHODS 

In Paper I we have developed a numerical code that can yield accurate solutions 
to the discrete form of the coagulation equation (equation [1]) with a reasonable 



number of numerical mass bins. A detailed description of the code can be found 
in Paper I. In this section we provide a brief summary of the algorithm and its 
numerical parameters. We also provide additional information on the accuracy of 
the numerical results. 

Our numerical code uses a combination of linearly and logarithmically spaced 
numerical mass bins. The first Nf,d numerical mass bins are linearly spaced with 
fhk = k. They have boundaries rhk±i/2 = k ± 1/2 and width Afhk = fhk+1/2 — 
™k-i/2 = 1- The next Nm x Ad ec numerical mass bins are logarithmically spaced, 
with Nbd bins per decade of mass; thus the kth mass is fhk = (fhk+1/2 + fhk-i/2)/2, 
with fh k+ i/ 2 /fh k _i/ 2 = lO 1 /^". There are in total A max = (A dec + 1)A M mass 
bins, and the mass of the most massive particles in the computational grid, m max , 
is approximately NbdlQ Ndcc • Initially, the number of "active" bins Ab; n <§C A max . 
At the end of each time step, Abi n is increased (if necessary) to include all bins with 
n mim where N k is the total number of particles in bin k and n m [ n is a numerical 
parameter; it is also increased if Ajv bin +i becomes comparable to the power-law 
extrapolation from Ajv bin . Before Abi n reaches A max , A max (or equivalently m max ) 
has no effects on the numerical results. The numerical parameter n m i n is specified 
in units of no and, e.g., n m i n = 10~ 30 in units of n is equivalent to n m i n = 1 and 
uq = 10 30 in physical units. Thus the effect of the numerical parameter n m i n is 
similar to not allowing fractionally occupied numerical mass bins to interact. 

The fundamental quantity evolved by our numerical code is the total mass M k 
in bin k. During a time step, the code calculates for each combination of i and 
j (with i < j < A b in) the mass loss from bins i and j due to mergers between 
particles in those bins and distributes the total mass of the merger products among 
the mass bins. Thus the code conserves mass exactly. For i < j < N(,d, it is correct 
to assume that the merging particles have masses fhi and fhj and that the merger 
products have mass fhi + fhj. For % < j and j > A^, we assume that the particles 
in bin i have mass fhi (which is exact for i < A^) and that the mass distribution 
within bin j follows a power-law distribution: 

Pj(m) = c j (m/fhj- 1 / 2 ) qj for mj-1/2 < m < fh j+1 / 2 (4) 

where pj(m)dm is the total mass of particles with mass between m and m+dm. The 
merger products have masses between fhi + fhj_i/ 2 and fhi + m J+ i/2, and they are 
either added to a single bin k (if fhk-1/2 < fhi + fhj_i/ 2 and fhi + fhj + i/ 2 < mfc+1/2) 
or distributed between bins k and k + 1. In equation (4), the power-law index qj is 
obtained from the masses in the adjacent bins, 
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and the normalization constant Cj from the constraint 

rrhj+i/2 

/ dm Pj(m) = Mj. 

Jrhj-l/2 



(6) 



Our numerical code uses the second-order Runge-Kutta method with a variable 
time step. The time step is continuously adjusted so that the fractional change of 
each Mk per time step is less than 5m and the mass loss from bin k does not exceed 



In Paper I we have compared in detail the numerical solutions from our code 
to the exact analytic solutions for Aij = 1, (i + j)/2, and ij, with the last case 
at t < t only. The accuracy of the numerical solutions is extremely insensitive to 
5m and n m ; n and improves rapidly with increasing Nm- Hereafter, unless stated 
otherwise, the numerical results were obtained using 5m = 5%, n m i n = 10~ 30 (in 
units of no), and Nm = 40. 

We report here several additional tests of our code. Ziff (1980) has constructed 
three forms of rate coefficients, with a parameter 7, for which a single moment 
■M 7 (t) = Y^kLi m 1 n k{i) of the mass spectrum rifc(t) can be calculated analytically. 
We have obtained numerical solutions for a few of these rate coefficients (including 
both orderly and runaway growth cases) and have confirmed that the numerical 
results for the moment M. y (t) agree with the analytic results, with accuracy similar 
to what was found for the three cases with exact analytic solutions. 

In Section 4.3 we shall be interested in extending the calculations for some of 
the runaway growth cases with v < 1 and A > 1 to t > to- Therefore, another 
test that we have performed is to extend the comparison for the case = ij to 
t > to- For = ij, the evolution of the mass spectrum at t > to depends on 
whether or not the runaway particle (gel) and the other particles (sol) interact. 
The (unmodified) coagulation equation is valid if there is no sol-gel interaction, and 
it has an exact analytic solution with nk(t) oc t -1 /c -5 / 2 at large k for all t > to 
(Leyvraz and Tschudi 1981, Ziff et al 1983). (The coagulation equation can be 
modified to take into account sol-gel interaction, and an exact analytic solution 
also exists for this modified coagulation equation; Ziff et al 1983.) Since our code 
does not take into account sol-gel interaction and does not allow merger products 
with masses greater than m max to interact, we expect the numerical solution at 
t > to to agree with the analytic solution to the unmodified coagulation equation, 
except for ~ m max . We have integrated the case A^ = ij up to t = 1.25 and 
have found that the numerical solution at t > to = 1 is in excellent agreement with 
the analytic solution for mk ^ 0.01m max - 

A quantity that will be discussed extensively in Section 4 is the logarithmic 
slope d In n/d In m of the mass spectrum. For the numerical results, we use 



d\nn I 'd\nm (mk) = qk — 1 



(7) 
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where qk is defined in equation (5). Equation (7) is consistent with the power- 
law approximation used by the code since under this approximation the number 
distribution of particles within the numerical mass bin k is n/-(m) = pk(m)/m oc 

To determine the accuracy of the numerical results for d In n/d In m, we compare 
the numerical and analytic results for the three cases with exact analytic solutions. 
The analytic solutions are of the form n oc m~ T exp[— b(t) m] or <ilnn/<ilnm = 
— t — b{t)m for m 3> 1, where r = 0, 3/2, and 5/2 for = 1, (i + j)/2, and 
ij, respectively. Thus n oc m~ T or <ilnn/<ilnm = — r for 1 <C m <C m* when the 
characteristic mass m*(£) defined in equation (9) is large. In figure 2(a) we show for 
each case the numerical and analytic d\nn/d\nm at m < m*(£) at a given time. As 
we noted in Paper I, there is a small lag in the evolution of the numerical solutions 
for Aij = (i + j)/2 and ij. Therefore, in these cases, the numerical results are 
compared to the analytic results at a slightly earlier time. There are fluctuations 
in the numerical results in the first decade of the logarithmically spaced mass bins 
(1 < m/N bd < 10) due to the discreteness of the mass bins, but the fluctuations 
are < 0.015. The numerical results are much smoother and much more accurate 
outside this mass range. We can see from figure 2(a) that r can be determined from 
the numerical results at 1 <C m <C m* to better than ±0.001. 

Figure 2(b) is similar to figure 2(a), but it shows the mass range m > m*(t). 
The differences between the analytic results, which decrease linearly with mass, 
and the numerical results are small, but there is a small curvature in the numerical 
results, and the numerical results become increasingly higher than the analytic 
results with increasing mass. (This is consistent with the observation in Paper I 
that the numerical solutions show a slightly slower exponential decay at the high- 
mass end of the mass spectrum.) As a result, if we fit the numerical results near 
d\nn/d\nm = —20 (or —30) to a straight line d In n/d In m = —9 — bm, the resulting 
values for 9 are greater than the correct values (which are r as given above) by 0.13- 
0.15 (or 0.29-0.36). This is the accuracy to which we can check whether a numerical 
solution is consistent with d In n/d In m = —9 — bm and a given ^ at m > m*. 

In Paper I we have discussed our numerical code in the context of numerical 
codes in the astrophysics literature. Other recent numerical codes include those by 
Krivitsky (1995), Hill and Ng (1996), and Tzivion et al (1999). As we mentioned 
in Section 1, the numerical code used by Krivitsky does not conserve mass and 
would have difficulty following the evolution of the mass spectrum accurately for 
a long time. The numerical code described by Hill and Ng conserves mass but 
uses a relatively simple algorithm for distributing merger products. We had in fact 
tried a similar algorithm for distributing merger products (Quinlan and Shapiro 
1989) before we developed the algorithm based on the power-law approximation 
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and had found that the high-mass end of the mass spectrum converges very slowly 
with increasing grid resolution (iV^) if the rate coefficient increases steeply with 
the mass of the particles (i.e., if v and/or A is large). Tzivion et al have developed a 
mass-conserving numerical code that evolves separately the total number (N k ) and 
mass (M k ) of particles in a numerical mass bin k. The numerical solutions obtained 
using this code appear to converge rapidly with increasing grid resolution for the 
case Am^m* oc m + m', but there was no demonstration that this is also true for rate 
coefficients with steeper mass dependence. 

4. RESULTS 

In this section we present the numerical solutions to the discrete form of 
the coagulation equation for the rate coefficient and initial conditions described 
in Section 2 (see also figure 1). For the cases with v < 1 (Sections 4.1-4.3), 
whenever possible, we compare the properties of the numerical solutions to the 
predictions from self-similar analysis. Hereafter, unless otherwise stated, the self- 
similar analysis results cited can be found in van Dongen and Ernst (1985a, 1988). 

4.1. Orderly Growth Cases with v < 1 and A < 1 

An example of the numerical results for the mass spectrum evolution for v < 1 
and A < 1 is shown in figure 3. 3 In this and all other cases with v < 1 and A < 1, 
the mass spectrum evolves in an orderly fashion. For these cases, we stopped the 
numerical integrations when the asymptotic behaviors of the solutions were clear 
and before iVbin reached iV max . 

For orderly growth with v < 1 and A < 1, self-similar analysis predicts that 
self-similar solutions have the form 

n k (t) = m*{t)~ 2 ip[m k /m*{t)} (8) 

where (p(x) is a scaling function, and the characteristic mass m*(t) scales as t 1 ^ 1- ^. 
Different definitions of m* (t) , which correspond to different scales for x = m k / (t) 
and (p(x), can be used. We adopt 

m*(t) = M 3 (t)/M 2 (t) (9) 

where M.t{t) = J^feLi m k n k{'t) is the £th moment of the mass spectrum. This choice 
of m* is convenient because it can also be used for runaway growth with v < 1 and 
A > 1 (Section 4.3). 

3 In figure 3 and all subsequent figures, the numerical mass spectrum plotted is 
n{rhk) = N k for k < N bd and n(m k ) = n k (fh k ) for k > N bd . 
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In all cases, the numerical solution tends toward a self-similar solution of the 
form (8) at large t. This is illustrated in figure 4, where we plot the numerical 
solution at three different times in the form of log m 2 n as a function of log m/m* 
for the case shown in figure 3. In the scaling limit (8), m 2 n = x 2 (p(x). We 
have evaluated the exponent z(ti) = ln[m*(£j+i)/m*(£j_i)]/ln(£j + i/£j_i) from the 
numerical results for m*(t) at output times U-i, ti, and ti+i for all cases with v < 1 
and A < 1. In all but two cases, the exponent z at large t agrees with 1/(1 — A) to 
better than one part in 1.5 x 10 3 . For the two cases with (fi, v) = (1/6,2/3) and 
(1/3, 1/2), the agreement between the numerical results at the end of the numerical 
runs, z = 5.985 and 5.963, and the analytic result, 1/(1 — A) = 6, is slightly worse, 
but the numerically determined exponents are still slowly increasing with time at the 
end of the numerical runs, indicating that the numerical results have not completely 
reached the asymptotic regime. 

Figure 5 shows the numerical results for the scaling function <p(x) for all cases 
with v < 1 and A < 1. The location of the peak of x 2 (p is not very sensitive to 
fi or z/, and it is at x = 0.33-0.71. Self-similar analysis predicts that the scaling 
function <p(x) decays exponentially at large x. For v < 1, <p(x) oc x~ e exp(— bx) 
or d In ip/d In x = —9 — bx, with 9 = A, at large x. The detailed behavior of <p(x) 
at large x for v = 1 depends on the specific form of A^. For A^ in equation 
(3) with v = 1, the large-x behavior of (p(x) is similar to that for v < 1, but 
9 = (/j + 3)/2 if -1 < n < (van Dongen 1987b; see also Ernst et al 1984). In 
all cases, the numerical <p(x) decays exponentially at large x, and d\mp / d\nx at 
large x is consistent with the analytic result (see Section 3 for a discussion of the 
accuracy of the numerical results at large x). 

The behavior of the scaling function (p(x) at small x is qualitatively different 
for \i < 0, jU = 0, and > 0. For the cases with \i < (figure 5(a)), (p(x) 
also decays exponentially at small x, because light particles are rapidly accreted 
by heavy particles. Self-similar analysis predicts that ip(x) oc x~ a exp(6x M //u) or 
d In ip/ d In x = —a + bx^ at small x, where o and b are constants that depend on 
the specific form of Aij. For A^ in equation (3), a = 2 if v > and a = 1 if v = 0. 
In figure 6, we plot the numerical results for dhup/dhi function of x^ for 

the cases with \i = —1/6 to show that the numerical results indeed have the form 
d\ii(f/d\nx = —a + bx^ at large x^ (or small x). Furthermore, in all cases with 
fx < 0, the value of a from least-squares fit is in agreement with the analytic result 
to better than ±0.004. 

For the cases with fx = 0, i.e., A^ = [i v + j")/2, the scaling function shows a 
power-law behavior ip(x) oc x~ T at small x (figure 5(b)). The numerical results for 
the exponent r are 0.000, 1.001, 1.033, 1.109, 1.216, 1.347, and 1.500 for v = 0, 1/6, 
1/3, 1/2, 2/3, 5/6, and 1, respectively. Self-similar analysis gives r = 2 —p\/w, 
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where p\ and w depend on the specific form of A^-. van Dongen and Ernst (1985b) 
have used this expression to derive analytic lower and upper bounds on r for Aij = 
[i v +j 1/ )/2. Our numerical results are consistent with these bounds. Note that the 
exponent r is discontinuous at v = 0: r decreases from 1.5 at v = 1 to 1 as v \ 0, 
but r = at v = (recall that the z/ = and z/ = 1 cases have exact analytic 
solutions). In contrast, for oc (i + j) u (which also has \i = 0), the exponent r 
decreases smoothly from 1.5 at v = 1 to at v = (van Dongen and Ernst 1985c, 
Krivitsky 1995). 

For the cases with \i > (figure 5(c)), the scaling function (p(x) oscillates 
around a power law at small x, with the fractional amplitude of the oscillation 
decreasing as x — > 0. This damped oscillatory approach to a power law is shown 
more clearly in figure 7, where we plot the logarithmic slope din ip /din x as a 
function of logx for the cases with \i = 1/6. In all cases, the logarithmic slope 
tends to a constant value as x — > 0, and the asymptotic value is consistent with 
the leading small-x behavior predicted by self-similar analysis: (p(x) oc a; _(1+A ) or 
dlmp/dlnx = — (1 + A). The oscillation appears to be periodic in the variable Inx, 
but it is difficult to determine this accurately because of the limited number of 
cycles seen in our numerical results. 

The damped oscillatory behavior at small x for \i > and A < 1 was not known 
previously. In their self-similar analysis, van Dongen and Ernst (1988) were unable 
to find higher-order corrections to the leading small- a; behavior of (p(x) for \x > and 
A < 1, and they raised the possibility that physically acceptable self-similar solutions 
may not exist. As we have just shown, there are indeed physically acceptable self- 
similar solutions, and they are reached from monodisperse initial conditions. Based 
on our numerical results, we suggest that the leading small- a; behavior and the first 
correction could be of the form (p(x) oc £ _(1+A )[1 + f(x) cos(£? Inx + C)], where 
f(x) is an increasing function of x, possibly Ax a with a > 0. The failure to find 
the first correction in the self-similar analysis is probably due to the unusual form 
of this correction. 

As we mentioned in Section 1, Krivitsky (1995) has obtained numerical 
solutions to the continuous form of the coagulation equation for A m:rn ' oc (mm') A / 2 
with A < 1. Krivitsky concluded that the numerical solutions are self-similar at large 
time but that the asymptotic behavior at small mass is not a power law. The latter 
conclusion is different from ours and, we believe, incorrect. We can understand 
why Krivitsky reached this conclusion by examining figure 5(b) of Krivitsky (1995), 
where the numerical results for d In n/d In m are shown for the case A = 0.4. By 
the last time shown, the numerical results have converged to a self-similar form for 
x > 10 -4 , and the logarithmic slope is indeed decreasing with decreasing x over 
the range 10 -4 < x < 10 _1 . However, as we can see from our figure 7, over this 
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range in a;, dlmp/dlnx is in fact decreasing from a maximum to a minimum in its 
oscillatory approach to a constant value. Therefore, the incorrect conclusion was 
reached because Krivitsky did not evolve the numerical solutions for a sufficiently 
long time to see the true asymptotic behavior at small mass. 

4.2. Orderly Growth Cases with v < 1 and A = 1 

On the borderline A = 1 and v < 1, the numerical mass spectrum evolves in 
an orderly fashion, but with the characteristic mass m* (t) increasing exponentially 
with time. For these cases, we set iVdec = 19 and stopped the numerical integrations 
as soon as iVbi n reached iV max , i.e., when the mass spectra extended over 20 decades 
in mass. 

We distinguish the cases with fx = and fx > 0. As we discussed above, the case 
fx = is Aij = (i + j)/2 with exact analytic solution, and the numerical solution 
for this case is in excellent agreement with the analytic solution (see Section 3 and 
figures 2 and 5(b)). The mass spectrum tends toward a self-similar solution of 
the form (8) at large t, with m*(t) oc e* (see figure 8), and the scaling function 
(f(x) oc x~ 3 / 2 at small x and oc x~ e exp(— bx) with 6 = (fi + 3)/2 = 3/2 at large x. 

For fx > 0, van Dongen and Ernst (1988) have derived a modified self-similar 
solution: 

"fcW = (™* lnm*)"V(m*;/m*) (10) 

where (lnm*) 2 = a + bt and a and b are constants. The scaling function <p(x) is 
predicted to scale as x~ 2 at small x and x~ x exp(— bx) at large x. 

The numerical results for m*(t) for the three cases with fx > (and also the 
case fx = 0) are shown in figure 8. In each case, we have fitted the numerical 
results at the last two output times to lnm* = a + bt and (lnm*) 2 = a + bt, and 
they are shown as dotted and solid lines, respectively. For (fx, u) = (1/2,1/2), 
(lnm*) 2 = a + bt provides a reasonably good fit to the numerical results at large 
t. For (fx, u) = (1/3,2/3) and, in particular, (1/6,5/6), the numerical results at 
large t show deviations from (lnm*) 2 = a + bt. The deviations in the last two 
cases are probably due to the numerical results not having completely reached 
the asymptotic regime, but we cannot rule out the possibility that the asymptotic 
behavior is different from the analytic prediction. 

If a numerical solution approaches the self-similar solution (10), we expect 
m 2 nlnm* — > x 2 ip(x) and d In n/d lnm — > dhup / 'dlnx. In figure 9 we show 
d In n/d lnm as a function of log m/m* at three different times for all cases with 
fx > and A = 1. For (fx, v) = (1/2, 1/2), d In n/d lnm has converged to d\mp/ 'dlnx 
at x = m/m* > 10 -5 , but the convergence at x < 10 -5 is very slow and is not 
complete by the end of the numerical run. The range of x over which d In n/d lnm 
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has converged by the end of the numerical run is wider for smaller n. (A similar 
analysis of m 2 nlnm t reveals a small increase in the normalization of m 2 nlnm, 
with time in the range of x where d In n/d In m has converged. This increase is 
more pronounced for smaller \i and is probably due to m* not having completely 
reached the asymptotic regime.) For (/i, v) = (1/6,5/6), it is reasonably clear that 
dlmp/dlnx — > —2 in the small- x limit, consistent with the analytic prediction. For 
(/U, v) = (1/3,2/3) and (1/2,1/2), the small-x behaviors are less certain because 
of the slow convergence at small x, but they also appear to be consistent with the 
analytic prediction. Finally, in all three cases, the large-x behavior of the numerical 
dlmp/ d\nx is consistent with the analytic prediction that d\mp/ d\nx = —1 — bx. 

4.3. Runaway Growth Cases with u < 1 and A > 1 

In all cases with v < 1 and A > 1, the numerical mass spectrum develops 
in a finite time to a power-law tail, nk oc /c _T , at large mass that violates mass 
conservation, and runaway growth is expected to start at t cr it = to i n the limit 
uq — > oo. For most cases with v < 1 and A > 1, we stopped the numerical 
integrations as soon as iVbin reached iV max and the mass spectra extended over 
20 decades in mass; so the numerical solutions approach very close to but do not 
exceed to- To study the transition at t = to, we extended the integrations for the 
cases (/U, z/) = (1/3,1) and (2/3,2/3) to t > t . Figure 10 shows the numerical 
results for the mass spectrum evolution for the case (,(/, v) = (1/3, 1) at t < t (solid 
lines) and t > to (dashed lines). 

For runaway growth with v < 1 and A > 1, self-similar analysis predicts that 
self-similar solutions close to but before to have the form 

n k (t) = m*(t)~ T Lp[m k /m*(t)] (11) 

where the scaling function <p(x) oc x~ T at small x, the characteristic mass m*{t) 
diverges as (to — t) -1 /' 7 , and a = A + 1 — r. With m*(t) diverging at to, the self- 
similar solution (11) has nfc(to) oc k~ T at large k. In all cases, the numerical solution 
tends toward a self-similar solution of the form (11) as t j to, and the numerical 
(p(x) is indeed a power law at small x. This is illustrated in figure 11, where we plot 
the numerical solution at three different times (< to) in the form of log(m 2 nm* -2 ) 
as a function of log m/m* for the case shown in figure 10, with the exponent r 
determined from the numerical solution itself (see table 1 and discussion below) . In 
the scaling limit (11), m 2 nml~ 2 = x 2 tp(x). 

Despite the change in the form of the self-similar solution, the analytic 
predictions for the large- a; behavior of (p(x) are similar to those for orderly growth 
in Sections 4.1-4.2: tp(x) oc x~ exp(— bx) or dlinp/dlnx = —9 — bx, where 9 = A 
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if v < 1, 9 = {n + 3)/2 if i/ = 1 and < \i < 1, and 9 = 5/2 if v = fi = 1. In all 
cases, the numerical (p(x) decays exponentially at large x, and d\n(p/d\nx at large 
x is consistent with the analytic result. 

The analytic predictions for the exponents r and a are r = (A + 3)/2 and 
a = X + 1 — r = (A — l)/2 (Hendriks et al 1983, van Dongen and Ernst 1985a, 
1988). For comparison, we have determined r and a from the numerical results 
for d In (p/d In x and m*(t), respectively. In figure 12 we show the numerical results 
for dliap/dlnx for the cases with A = 4/3. The numerical results clearly converge 
to constant values at small x. The constant asymptotic values are consistent with 
<p(x) oc x~ T or d\mp/ d\nx = —t at small x and directly yield the values of r. It is 
also clear from figure 12 that the asymptotic values and hence r for these cases with 
the same A(= 4/3) are different from each other and from the analytic prediction 
that d In ip/ d In x = — r = — (A + 3)/2 = —13/6 (dashed line in figure 12). The 
numerical results for the exponent r for all cases with v < 1 and A > 1 are listed 
in table 1. The only case where the exponent r agrees with the analytic prediction 
(A + 3)/2 is the case v = ji = 1 (i.e., the case Aij = ij with exact analytic solution). 
In all other cases, the exponent r is less than (A + 3)/2. For a given A, the deviation 
of t from (A + 3)/2 is largest for v = 1 and smallest for v = ji. 

To determine the exponent a, we fit the numerical results for m*(t) at output 
times ti, and t i+ i to m*(i) = C(t' — t) -1 / " to obtain C(tj), t' (ti), and a'(ti). 
In most cases, cr'(t) has converged to a constant value by the end of the numerical 
run and directly yields a. In the remaining cases, we extrapolated cr'(t) to the 
limit m*(t) — > oo to obtain a, but the difference between cr'(t) at the end of the 
numerical run and a is < 0.002. The numerical results for the exponent a are listed 
in table 1, together with A + 1 — r. Except for the case v = ji = 1, the exponent 
a is greater than the analytic prediction (A — l)/2. For a given A, the deviation 
of a from (A — l)/2 is largest for v = 1 and smallest for v = ji. We note that the 
numerical results for a and r are consistent with one another in that they satisfy 
the relation a = A + l — r for self-similar solutions of the form (11) to within ±0.001. 

The procedure described above for determining the exponent a also yields the 
time to- The numerical results for to, in units of l/(no^4n), are listed in table 1. 
Since we expect to ~ 1/[(A — l)no^4n] (see Paper I), it is convenient to parameterize 
t as t = — l)n in]. The parameter K is shown in figure 13. We find that 

K = 1-2 if A > 1.1 and that, for a given A, K is smallest for v = 1 and largest for 
v = fx. For v = 1, the parameter K shows a maximum at A ~ 1.3. For v = fj,, K 
increases monotonically with decreasing A, but it is unclear whether K approaches a 
constant value or diverges as A — > 1. Finally, we note that the numerical results for 
to are consistent with the bound to > 1/[(A — l)no^4n] for in equation (3) and 
monodisperse initial conditions and the stronger bound to > 1/[(2 A_1 — l)no^4n] 
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for v = fj,, derived analytically by Hendriks et al (1983). 

We have found that the exponents r and a are in general different from the 
analytic predictions. Since the analytic prediction for a follows from that for r 
and the relation a = A + 1 — r for self-similar solutions of the form (11), and the 
numerical results for a and r satisfy this relation, we have essentially a discrepancy 
in the exponent r. Let us examine the derivation of the analytic prediction for r, 
which can be summarized as follows (see Hendriks et al 1983, van Dongen and Ernst 
1985a, 1988 for details). The mass flux from particles of mass mi < rrik to particles 
of mass rrii > nik is 

k oo 

F k(t) = Yl Yl rniAijmWrijit). (12) 

i=l j=k+l— i 

It is assumed that solutions to the coagulation equation at t > t violate mass 
conservation by having a non-zero and finite mass flux to the infinite-mass bin (i.e., 
in the limit k — > oo), which is possible only if the mass spectrum at t > to is of 
the form n^(t) oc k~ T at large k. With rifc(t) oc k~ T at large k, the mass flux 
F^it) oc k x+3 ~ 2r at large k. Thus r' = (A + 3)/2 if the mass flux -Ffc(t) is required 
to be non-zero and finite in the limit k — > oo. Since the self-similar solution (11) has 
nk(to) oc k~ T at large k, r = r' = (A + 3)/2 if we assume that the large-/c behavior 
at t > t is also valid at t = to. It should be noted that the arguments leading 
to t = (A + 3)/2 are not rigorous. In particular, as van Dongen and Ernst (1988) 
pointed out, one cannot exclude the possibility that the mass flux diverges at one 
instant of time, i.e., to- We have found that r < (A + 3)/2 in general. This implies 
that the mass flux at to? which is oc k x+3 ~ 2r at large /c, diverges as k — >■ oo. 
Thus the numerical solutions violate mass conservation at to by having a diverging 
mass flux to the infinite-mass bin. 

Since the mass flux cannot diverge for all times t > to, do the solutions at 
t > t have the form nk(t) oc &;~( A + 3 )/ 2 at large k for non-zero and finite mass flux 
to the infinite-mass bin? If so, how can this large-/c behavior at t > to be reconciled 
with that at t = to? To answer these questions, we have extended the numerical 
integrations for the cases (/U, v) = (1/3, 1) and (2/3,2/3) to t > to- The results for 
the mass spectrum evolution at t > t for the case (fx, v) = (1/3, 1) are shown as 
dashed lines in figure 10. As in the test case Aij = ij discussed in Section 3, the 
numerical solution at mk ~ m max is affected by the finite value of the maximum 
particle mass, m max , in the computational grid. Otherwise, the mass spectrum at 
t > t is indeed of the form n^t) oc /c~( A + 3 )/ 2 at large k, with the value of the 
exponent confirmed by an analysis of the logarithmic slope. Note, however, that 
the range of k where nk(t) oc &;~( A + 3 )/ 2 shrinks toward infinity as t j to- The 
numerical solution for the case (fx, v) = (2/3, 2/3) shows the same large-/c behaviors 
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at t > to. Therefore, we conclude that the transition at t = to is accomplished as 
follows. As t f to, the solution tends toward a self-similar solution of the form (11), 
with rik(t) oc k~ T for 1 <C <C to*(£), r < (A + 3)/2 in general, and m*(£) — > oo. 
As t 1 1 , nfe(t) oc /c"( A + 3 )/ 2 for m fc > m*(t) and m'*(t) -> oo. 

4.4. Runaway Growth Cases with v > 1 

Examples of the numerical results for the mass spectrum evolution with n m i n = 
10 -30 (solid lines) and 10 -35 (dotted lines) for v > 1 are shown in the lower panels 
of figure 14. In all cases with v = 2 and 3/2 and the three cases with v = 7/6 
and /i < (see, e.g., figure 14(a)), the behaviors of the numerical solutions for both 
^min = 10~ 30 and 10~ 35 are qualitatively similar to those found in Paper I for other 
forms of Aij with v > 1. In the early stages, the mass spectrum appears to decay 
exponentially at large mass. After some evolution, the mass spectrum shows an 
exponential drop in an intermediate mass range and a power-law tail of the form 
tik oc k~ v (or n oc m~ u ) at large mass. For the three cases with v = 7/6 and fx > 1/3, 
the numerical solutions with n m i n = 10 -30 do not develop the m~ v tail when iVb; n 
reaches A max (see, e.g., figure 14(b)), and their behaviors are qualitatively similar 
to those for the runaway growth cases with v < 1 and A > 1 (Section 4.3; figure 
10). When n m i n is decreased to 10~ 35 , the \i = 1/3 case does develop an m~ u tail 
when A^ bin reaches A max (figure 14(b)), but the other two cases do not. It is not 
feasible to compute solutions for much smaller n m i n , but we strongly suspect that 
the remaining two cases would develop the m~ v tail for sufficiently small n m ; n . 

The fact that the tail at large mass is of the form n oc m~ v or dlnn/dlnm = —v 
is illustrated in the upper panels of figure 14, where we plot the numerical results 
for d In n/d In m at the specified times, just after the formation of the tail for the 
runs with n m i n = 10 -35 . As found in Paper I, the time at which the power-law 
tail develops decreases slowly toward zero as n m ; n decreases (lower panels of figure 
14). We do not repeat the demonstrations here, but it was shown in Paper I that 
numerical solutions with different maximum particle mass, m max , included in the 
computational grid are identical in the overlapping mass range and that the power- 
law tail simply extends to larger particle mass m when m max is increased. (It was 
also shown that the numerical solutions converge by = 40 and 5m = 5%.) 
Therefore, in the limit n m ; n — > and m max — > oo, the numerical solutions for all 
cases with v > 1 (with the possible exception of the cases with v = 7/6 and ji > 1/3) 
should develop in an infinitesimal time power-law tails of the form n oc m~ v that 
extend to arbitrarily large m max . However, this mass spectrum is not self-consistent 
because the power-law tail implies a mass flux and, if 1 < v < 2, a cumulative mass 
that diverge with m max (see footnote 2 and Paper I). Thus, as in Paper I, the 
numerical results strongly suggest that there are no self-consistent solutions to the 
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coagulation equation at any time if v > 1. Furthermore, since the formation of 
the power-law tail in the coagulation equation solution with finite n m i n probably 
corresponds to the onset of runaway growth in Monte Carlo simulations with finite 
no (see Paper I) and decreasing n m i n is equivalent to increasing no (see Section 3), 
the time t C rit for the onset of runaway growth for v > 1 should decrease slowly 
toward zero with increasing no- 

5. CONCLUSIONS 

We have conducted a systematic survey of numerical solutions to the 
coagulation equation (1) for a rate coefficient of the form Aij oc {i^j v + i v j^) 
and monodisperse initial conditions. The numerical results confirmed that there 
are three classes of rate coefficients with qualitatively different solutions to the 
coagulation equation. 

For v < 1 and A = \x-\-v < 1 (Section 4.1), the numerical solution evolves in an 
orderly fashion and tends toward a self-similar solution of the form (8) at large time 
t. The time dependence of the characteristic mass m*(t) at large t and the behaviors 
of the scaling function <p(x) at small and large x = nik/m*(t) agree with the analytic 
predictions of van Dongen and Ernst (1985a, 1988). In particular, for the subset of 
cases with \i > 0, we disagreed with the earlier numerical study by Krivitsky (1995) 
and found that <p(x) does in fact approach the analytically predicted power-law 
behavior <p(x) oc x~( 1+x ^ at small x, but in a damped oscillatory fashion that was 
not known previously (figure 7). For fi = 0, we determined the exponent r of the 
power-law behavior f(x) oc x~ T at small x. 

On the borderline v < 1 and A = 1 (Section 4.2), the numerical solution 
evolves in an orderly fashion. For fx = 0, i.e., oc i + j, the numerical solution 
is in excellent agreement with the exact analytic solution and tends toward a self- 
similar solution of the form (8) at large t. For a > 0, the numerical solution appears 
to tend toward a self-similar solution of the form (10), but we had limited success 
in comparing the behaviors of m*(t) and <p(x) at small x to the analytic predictions 
because the convergence to the self-similar solution is very slow. 

For v < 1 and A > 1 (Section 4.3), the numerical mass spectrum develops 
in a finite time to a power-law tail, oc /c _T , at large k that violates mass 
conservation, and runaway growth/gelation is expected to start at t CT i t = to in 
the limit that the initial number of particles no — > oo. As t | to, the numerical 
solution tends toward a self-similar solution of the form (11), with (p(x) oc x~ T 
at small x and m*{i) diverging as (to — t) -1 ^ . The exponent r is in general 
less than the analytic prediction (A + 3)/2, and the exponent a greater than the 
analytic prediction (A — l)/2, but they satisfy the relation a = A + 1 — r for self- 
similar solutions of the form (11) (table 1). We studied the dependence of to on the 
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exponents fj,, v, and A and found that to = K/[(\ — l)noAu] and K = 1-2 if A > 1.1 
(table 1; figure 13). At t > t , n k oc &;-( A + 3 )/ 2 for m k > m'*(t), with m^(t) — > oo 
as t i t . 

For v > 1 (Section 4.4), the behaviors of the numerical solution are qualitatively 
similar to those found in Paper I: the numerical mass spectrum develops a power-law 
tail of the form n k oc k~ v at large k that is not self-consistent, and the time at which 
the power-law tail develops decreases toward zero as the numerical parameter n m i n 
decreases. The numerical results strongly suggest that there are no self-consistent 
solutions to the coagulation equation at any time and that runaway growth/gelation 
is instantaneous in the limit no — > oo. They also indicate that the time t cr it, in units 
of l/(?ioAii), for the onset of runaway growth decreases slowly toward zero with 
increasing no, consistent with recent Monte Carlo simulation results (Malyshkin 
and Goodman 2001). 

The results presented in this paper and Paper I suggest several problems 
for future investigations. First, as we pointed out in Section 4.1, the failure to 
find the first correction to the leading small- x behavior (p(x) oc a;~( 1+A ) for the 
orderly growth cases with fx > and A < 1 in previous self-similar analysis is 
probably due to the unusual, damped oscillatory form of this correction. Given the 
information provided by the numerical results, it may now be possible to derive 
the first correction analytically. We suggested that the first correction could, e.g., 
be of the form x~( 1+x ^ f(x) cos(Blnx + C), where f(x) is an increasing function 
of x. Second, we have found that the exponent r for the runaway growth cases 
with v < 1 and A > 1 is in general different from existing analytic prediction. 
It is important to investigate whether analytic prediction matching the numerical 
results could be derived. In this connection, it would be useful to obtain numerical 
solutions for other forms of Aij and determine whether the exponent depends on 
the specific form of A^. Finally, it should be emphasized that rate coefficients 
with v > 1 do arise, and are of great interest, in astrophysics (see, e.g., Lee 1993, 
2000, and references therein). For astrophysical applications, we are interested in 
systems with finite no and interactions with the runaway particle. Thus, for v > 1, 
a detailed comparison of numerical solutions with finite n m ; n and Monte Carlo 
simulations with finite no should be conducted to test the correspondence between 
them, and the question of whether the coagulation equation can be modified to take 
into account the interactions between the runaway particle and the other particles 
should also be investigated. 
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Table 1. 


Runaway growth 


cases with v 


< 1 and A > 1 




A 




V 


T 


a 


A + 1 - T 


Tin Aa i tn 

' U 11 u 


7/6 


1/6 


1 


2.012 


0.154 


0.155 


7.136 


7/6 


1/3 


5/6 


2.038 


0.128 


0.129 


9.146 


7/6 


1/2 


2/3 
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2.054 


0.112 


0.113 


10.633 
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7/12 


7/12 
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0.110 


0.110 


10.855 


4/3 


1/3 
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2.076 


0.257 


0.257 


3.741 


4/3 


1/2 


5/6 


2.103 


0.230 


0.230 


4.284 


4/3 


2/3 


2/3 


2.112 


0.221 


0.221 


4.492 


3/2 


1/2 


1 


2.166 


0.334 


0.334 


2.451 


3/2 


2/3 


5/6 


2.184 


0.316 


0.316 


2.639 


3/2 


3/4 


d/4 


Z. loD 


0.314 


0.314 


2.664 


5/3 


2/3 


1 


2.269 


0.398 


0.398 


1.750 


5/3 


5/6 


5/6 


2.276 


0.390 


0.391 


1.808 


11/6 


5/6 


1 


2.380 


0.453 


0.453 


1.307 


11/6 


11/12 


11/12 


2.381 


0.451 


0.452 


1.317 


2 


1 


1 


2.500 


0.500 


0.500 


1.000 
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Figure 1. The exponents \i and v of the cases for which we have obtained 
numerical solutions to the discrete form of the coagulation equation with rate 
coefficient Aij = (i^j" + i u j^)/2. The orderly growth cases with v < 1 and A < 1 
and on the borderline v < 1 and A = 1 are discussed in Sections 4.1 and 4.2, 
respectively. The runaway growth cases with v < 1 and A > 1 are discussed in 
Section 4.3, and those with v > 1 are discussed in Section 4.4. 
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Figure 2. (a) Comparison of numerical (solid lines) and analytic (dotted 
lines) results for d In n/d In m at m < m*(t). The lines, from top to bottom, are 
<ilnn/<ilnm + C for A^- = 1 at t = 10 8 , Ay = (i + j)/2 at t = 18, and = ij at 
t = 1; the constant C = 0, 1.45, and 2.4, respectively. For = (i + j)/2 and ij, 
there is a small lag in the evolution of the numerical solutions, and the numerical 
results are compared to the analytic results at a slightly earlier time (1 — e)t, where 
e = 3.5 x 10~ 4 and 1.1 x 10~ 4 , respectively, (b) Same as (a), but at m > m*(t). The 
lines, in decreasing steepness, are d In n/d In m for Aij = 1 at t = 10 8 , Aij = (i+j)/2 
at t = 18, and = ij at t = 1. 
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Figure 3. Evolution of the mass spectrum for (fj,, u) = (1/3,1/3). We plot 
log m 2 n as a function of log m since m 2 n is the total mass per unit logarithmic mass 
interval: J m 2 ndlnm = J mndm. 
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Figure 4. Approach to self-similar solution as t — > oo for (fj,, i>) = (1/3, 1/3). 
The numerical solution is plotted for t = 10 3 (solid line), 10 4 (dashed line), and 10 5 
(dotted line) in the form of log m 2 n as a function of log m/m*. 
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Figure 5. (a) Scaling function cp(x) for the cases with \x < and v < 1. 
The solid (dotted) lines in increasing width are for fi = —1/2 (—1/6) and v = 
0,1/3,2/3,1. (b) (p(x) for the cases with \x = and v < 1. The solid lines in 
increasing width are for v = 0, 1/6, . . . , 1. (c) (/'(x) for the cases with > and 
A < 1. The solid lines in increasing width are for = 1/6 and v = 1/6, 1/3, 1/2, 2/3. 
The dotted lines, offset vertically by —10, are for fx = 1/3 and v = 1/3 and 1/2. 
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Figure 6. Logarithmic slope din (p/ din x as a function of for the cases with 
\i = —1/6 and v = 0, 1/3, 2/3, 1 (solid lines from top to bottom). The dotted lines 
are the asymptotes approached by the numerical results at large x M (or small x). 
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Figure 7. Logarithmic slope dlmp/dlnx for the cases with \i = 1/6 and 
v = 1/6, 1/3, 1/2,2/3 (solid lines from top to bottom). The dashed lines show the 
leading small- a; behavior predicted by self-similar analysis: dlmp/dlnx = — (1 + A). 
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Figure 8. Characteristic mass m*(£) for the cases with A = 1 and v < 1. The 
points, from left to right, are the numerical results for \x = 0, 1/6, 1/3, and 1/2. The 
dotted and solid lines show lnm* = a + bt and (lnm*) 2 = a + bt fitted to numerical 
results at the last two output times. 
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Figure 9. Logarithmic slope dlnn/dlnm for the cases with \x > and A = 1. 
The top lines are the numerical results for (//, v) = (1/6, 5/6) at t = 47 (solid line), 
54 (dashed line), and 61 (dotted line), offset vertically by C = 1. The middle lines 
are for (fx, v) = (1/3, 2/3) at t = 63, 74, and 87, with C = 0.5, and the bottom lines 
are for (ji, v) = (1/2, 1/2) at t = 71, 86, and 102, with C = 0. 
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Figure 10. Evolution of the mass spectrum for (fj,, u) = (1/3,1) at t < to 
(solid lines) and t > t (dashed lines). 
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Figure 11. Approach to self-similar solution as t | to for (/z, i/) = (1/3,1). 
The numerical solution is plotted for t = 3.73 (solid line), 3.74 (dashed line), and 
3.7409 (dotted line) in the form of log(m 2 nm£~ 2 ) as a function of log m/m*. 
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Figure 12. Logarithmic slope dlnip/dlnx for the cases with A = 4/3. The 
solid lines from top to bottom are the numerical results for is) = (1/3,1), 
(1/2,5/6), and (2/3,2/3), and the dashed line is the analytic prediction for the 
small- a; behavior: dlmp/dlnx = — t = — (A + 3)/2 = —13/6. 




Figure 13. Parameter K = (A — l)no^4n£o for the cases with v < 1 and A > 1 
listed in table 1. The cases with v = 1 and v = [i are represented by squares and 
triangles, respectively, and the remaining cases are represented by crosses. 
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Figure 14. Evolution of the mass spectrum (lower panels) and logarithmic 
slope d In n/d In m at the specified time (upper panels) for |it = 1/3 and v equal 
to (a) 2 and (b) 7/6. The numerical solutions with n min = 10~ 30 (solid lines) 
and 10~ 35 (dotted lines) are shown. In the upper panels, the dashed lines indicate 
dlnn/dliim = —v. 



